knitr::opts_chunk$set(echo = TRUE)

cmdstanr::set_cmdstan_path(path = "C:/Users/kueng/.cmdstan/cmdstan-2.35.0")

library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
library(see)
library(beepr)
library(DHARMa)
library(digest)



source(file.path('Functions', 'ReportModels.R'))
source(file.path('Functions', 'PrettyTables.R'))
source(file.path('Functions', 'ReportMeasures.R'))
source(file.path('Functions', 'PrepareData.R'))

report_function_hash <- digest::digest(summarize_brms)
system("shutdown /a")
## [1] 1116
# Set options for analysis
use_mi = FALSE
shutdown = FALSE
report_ordinal = FALSE
report_hurdle = TRUE
do_priorsense = FALSE

options(
  dplyr.print_max = 100, 
  brms.backend = 'cmdstan',
  brms.file_refit = ifelse(use_mi, 'never', 'on_change'),
  brms.file_refit = 'on_change',
  #brms.file_refit = 'always',
  error = function() {
    beepr::beep(sound = 5)
    if (shutdown) {
      system("shutdown /s /t 180")
      quit(save = "no", status = 1)
    }
  }
  , es.use_symbols = TRUE
)


####################### Model parameters #######################

iterations = 2000
warmup = 1000
corstr = 'ar'
#corstr = 'cosy_couple'
#corstr = 'cosy_couple:user'


################################################################

if (corstr == 'ar') {
  autocor_str <- ~ ar(time = day, gr = coupleID:userID, p = 1)
  autocor_prior <- brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1) 
  suffix = paste0('ar_', as.character(iterations))
} else if (corstr == 'cosy_couple') {
  autocor_str <- ~ cosy(gr = coupleID)
  autocor_prior <- brms::set_prior("cauchy(0, 5)", class = "cosy", lb = -1, ub = 1)
  suffix = paste0('cosy_couple_', as.character(iterations))
} else if (corstr == 'cosy_couple:user') {
  autocor_str <- ~ cosy(time = day, gr = coupleID:userID)
  autocor_prior <- brms::set_prior("cauchy(0, 5)", class = "cosy", lb = -1, ub = 1)
  suffix = paste0('cosy_couple:user_', as.character(iterations))
}
df <- openxlsx::read.xlsx(file.path('long.xlsx'))
df_original <- df

df_double <- prepare_data(df, recode_pushing = TRUE, use_mi = use_mi)[[1]]

Constructing scales Re-coding pusing reshaping data (4field) centering data within and between

Modelling

Modelling

# For indistinguishable Dyads
model_rows_fixed <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb'
  )


model_rows_fixed_ordinal <- c(
  model_rows_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed[2:length(model_rows_fixed)]
)

model_rows_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]',
  'sigma'
)

model_rows_random_ordinal <- c(model_rows_random,'disc')
# For indistinguishable Dyads
model_rownames_fixed <- c(
    "Intercept", 
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily persuasion experienced", 
    "Daily persuasion utilized (partner's view)", # OR partner received
    "Daily pressure experienced", 
    "Daily pressure utilized (partner's view)", 
    "Daily pushing experienced", 
    "Daily pushing utilized (partner's view)", 
    "Day", 
    "Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean persuasion experienced", 
    "Mean persuasion utilized (partner's view)", 
    "Mean pressure experienced", 
    "Mean pressure utilized (partner's view)", 
    "Mean pushing experienced", 
    "Mean pushing utilized (partner's view)", 
    "Mean weartime"
  )


model_rownames_fixed_ordinal <- c(
  model_rownames_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed[2:length(model_rownames_fixed)]
)

model_rownames_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  "sd(Daily persuasion experienced)", 
  "sd(Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Daily pressure experienced)", 
  "sd(Daily pressure utilized (partner's view))", 
  "sd(Daily pushing experienced)", 
  "sd(Daily pushing utilized (partner's view))", 
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]',  
  'sigma'
)

model_rownames_random_ordinal <- c(model_rownames_random,'disc')
rows_to_pack <- list(
  "Within-Person Effects" = c(2,9),
  "Between-Person Effects" = c(10,16),
  "Random Effects" = c(17, 23), 
  "Additional Parameters" = c(24,25)
  )


rows_to_pack_ordinal <- list(
  "Intercepts" = c(1,6),
  "Within-Person Effects" = c(2+5,9+5),
  "Between-Person Effects" = c(10+5,16+5),
  "Random Effects" = c(17+5, 23+5), 
  "Additional Parameters" = c(24+5,25+6)
  )

HURDLE MODELS

# For indistinguishable Dyads
model_rows_fixed_hu <- c(
    'Intercept', 
    'hu_Intercept',
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb',
    
    # HURDLE MODEL
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'hu_persuasion_self_cw', 
    'hu_persuasion_partner_cw', 
    'hu_pressure_self_cw', 
    'hu_pressure_partner_cw', 
    'hu_pushing_self_cw', 
    'hu_pushing_partner_cw', 
    'hu_day', 
    'hu_weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'hu_persuasion_self_cb',
    'hu_persuasion_partner_cb',
    'hu_pressure_self_cb',
    'hu_pressure_partner_cb',
    'hu_pushing_self_cb',
    'hu_pushing_partner_cb',
    'hu_weartime_self_cb'
  )

model_rows_fixed_hu_ordinal <- c(
  model_rows_fixed_hu[1:2],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed_hu[3:length(model_rows_fixed_hu)]
)


model_rows_random_hu <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(hu_Intercept)',
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # HURDLE
  'sd(hu_persuasion_self_cw)',
  'sd(hu_persuasion_partner_cw)',
  'sd(hu_pressure_self_cw)',
  'sd(hu_pressure_partner_cw)',
  'sd(hu_pushing_self_cw)',
  'sd(hu_pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]',
  'sigma'
)

model_rows_random_hu_ordinal <- c(model_rows_random_hu,'disc')
# For indistinguishable Dyads
model_rownames_fixed_hu <- c(
    "Intercept", 
    "Hurdle Intercept",
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily persuasion experienced", 
    "Daily persuasion utilized (partner's view)", # OR partner received
    "Daily pressure experienced", 
    "Daily pressure utilized (partner's view)", 
    "Daily pushing experienced", 
    "Daily pushing utilized (partner's view)", 
    "Day", 
    "Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean persuasion experienced", 
    "Mean persuasion utilized (partner's view)", 
    "Mean pressure experienced", 
    "Mean pressure utilized (partner's view)", 
    "Mean pushing experienced", 
    "Mean pushing utilized (partner's view)", 
    "Mean weartime",
    
    # HURDLE
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Hu Daily persuasion experienced", 
    "Hu Daily persuasion utilized (partner's view)", # OR partner received
    "Hu Daily pressure experienced", 
    "Hu Daily pressure utilized (partner's view)", 
    "Hu Daily pushing experienced", 
    "Hu Daily pushing utilized (partner's view)", 
    "Hu Day", 
    "Hu Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Hu Mean persuasion experienced", 
    "Hu Mean persuasion utilized (partner's view)", 
    "Hu Mean pressure experienced", 
    "Hu Mean pressure utilized (partner's view)", 
    "Hu Mean pushing experienced", 
    "Hu Mean pushing utilized (partner's view)", 
    "Hu Mean weartime"
  )



model_rownames_fixed_hu_ordinal <- c(
  model_rownames_fixed_hu[1:2],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed_hu[3:length(model_rownames_fixed_hu)]
)



model_rownames_random_hu <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(Hurdle Intercept)', 
  "sd(Daily persuasion experienced)", 
  "sd(Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Daily pressure experienced)", 
  "sd(Daily pressure utilized (partner's view))", 
  "sd(Daily pushing experienced)", 
  "sd(Daily pushing utilized (partner's view))", 
  
  # Hurdle
  "sd(Hu Daily persuasion experienced)", 
  "sd(Hu Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Hu Daily pressure experienced)", 
  "sd(Hu Daily pressure utilized (partner's view))", 
  "sd(Hu Daily pushing experienced)", 
  "sd(Hu Daily pushing utilized (partner's view))", 
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]',
  'sigma'
)

model_rownames_random_hu_ordinal <- c(model_rownames_random_hu,'disc')
rows_to_pack_hu <- list(
  "Conditional Within-Person Effects" = c(3,10),
  "Conditional Between-Person Effects" = c(11,17),
  
  "Hurdle Within-Person Effects" = c(18,25),
  "Hurdle Between-Person Effects" = c(26,32),
  
  "Random Effects" = c(33, 46), 
  "Additional Parameters" = c(47,48)
  )

rows_to_pack_hu_ordinal <- list(
  "Intercepts" = c(1,7),
  "Conditional Within-Person Effects" = c(3+5,10+5),
  "Conditional Between-Person Effects" = c(11+5,17+5),
  
  "Hurdle Within-Person Effects" = c(18+5,25+5),
  "Hurdle Between-Person Effects" = c(26+5,32+5),
  
  "Random Effects" = c(33+5, 46+5), 
  "Additional Parameters" = c(47+5,48+6)
  )

Self-Reported MVPA

range(df_double$pa_sub, na.rm = T) 
## [1]   0 720
hist(df_double$pa_sub, breaks = 40) 

hist(log(df_double$pa_sub+00000000001), breaks = 40)

Hurdle Lognormal Model

formula <- bf(
  pa_sub ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID),
  
  hu = ~ persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  , autocor = autocor_str
) 

prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
  , brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA
  , brms::set_prior("normal(0.5, 2.5)", class = "Intercept", dpar = 'hu') # hurdle part
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
  , autocor_prior
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = hurdle_lognormal()
#)

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::hurdle_lognormal(), 
  #family = brms::hurdle_negbinomial(), 
  #family = brms::hurdle_poisson(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 42,
  file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
pa_sub_digest <- digest::digest(pa_sub)
check_brms(pa_sub, log_pp_check = TRUE)
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated possible pathological behavior:
##   Chain 1: E-BFMI = 0.011
##   Chain 2: E-BFMI = 0.102
##   Chain 3: E-BFMI = 0.035
##   Chain 4: E-BFMI = 0.065
## E-BFMI below 0.2 indicates you may need to reparameterize your model.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 1655 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## 
## Computed from 4000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -13342.0 237.3
## p_loo      4248.5 118.6
## looic     26684.1 474.6
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.0, 1.9]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     2081  55.7%   4       
##    (0.7, 1]   (bad)       439  11.8%   <NA>    
##    (1, Inf)   (very bad) 1216  32.5%   <NA>    
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(pa_sub, integer = TRUE, outliers_type = 'bootstrap')
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 77, observations = 3736, p-value < 2.2e-16
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000267666 0.003211991
## sample estimates:
## outlier frequency (expected: 0.00158993576017131 ) 
##                                         0.02061028
# summarize with rope range for hurdle part
summary_pa_sub_hurdle <- summarize_brms(
  pa_sub, 
  stats_to_report = c('pd','ROPE'),
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T) 
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
## Warning in summarize_brms(pa_sub, stats_to_report = c("pd", "ROPE"), rope_range
## = c(-0.18, : Coefficients were exponentiated. Double check if this was
## intended.
# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub$data$pa_sub[pa_sub$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_sub_continuous <- summarize_brms(
  pa_sub, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  rope_range = rope_range_continuous,
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T) 
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
## Warning in summarize_brms(pa_sub, stats_to_report = c("CI", "SE", "pd", :
## Coefficients were exponentiated. Double check if this was intended.
# Replace only the ROPE and % in Rope columns for rows with 'Hu'
summary_pa_sub <- summary_pa_sub_continuous

columns_to_replace <- c("ROPE", "inside ROPE")

summary_pa_sub[grepl('Hu', rownames(summary_pa_sub)), columns_to_replace] <- 
  summary_pa_sub_hurdle[grepl('Hu', rownames(summary_pa_sub_hurdle)), columns_to_replace]

# Print the updated dataframe
summary_pa_sub %>%
  print_df(rows_to_pack = rows_to_pack_hu)
exp(Est.) SE 95% CI pd ROPE inside ROPE Rhat Bulk_ESS Tail_ESS
Intercept 44.10*** 2.74 [38.98, 49.78] 1.000 [0.92, 1.08] 0.000 1.003 1322 2098
Hurdle Intercept 0.89 0.26 [ 0.51, 1.67] 0.654 [0.84, 1.20] 0.433 1.016 269 746
Conditional Within-Person Effects
Daily persuasion experienced 0.97 0.02 [ 0.93, 1.01] 0.938 [0.92, 1.08] 0.995 1.001 2240 2914
Daily persuasion utilized (partner’s view) 1.00 0.02 [ 0.96, 1.04] 0.504 [0.92, 1.08] 1.000 1.001 2365 2431
Daily pressure experienced 1.06 0.05 [ 0.96, 1.16] 0.882 [0.92, 1.08] 0.671 1.000 2068 2414
Daily pressure utilized (partner’s view) 0.99 0.05 [ 0.90, 1.08] 0.615 [0.92, 1.08] 0.878 1.000 2121 2895
Daily pushing experienced 1.03 0.03 [ 0.97, 1.09] 0.845 [0.92, 1.08] 0.942 1.002 2169 2791
Daily pushing utilized (partner’s view) 0.99 0.03 [ 0.93, 1.05] 0.657 [0.92, 1.08] 0.977 1.002 2410 2523
Day 1.13 0.08 [ 1.00, 1.29] 0.972 [0.92, 1.08] 0.240 1.001 2228 2579
Daily weartime NA NA NA NA NA NA NA NA NA
Conditional Between-Person Effects
Mean persuasion experienced 1.07 0.11 [ 0.87, 1.31] 0.746 [0.92, 1.08] 0.472 1.001 1442 1926
Mean persuasion utilized (partner’s view) 1.14 0.11 [ 0.94, 1.40] 0.908 [0.92, 1.08] 0.282 1.000 1299 1813
Mean pressure experienced 1.39** 0.17 [ 1.10, 1.79] 0.997 [0.92, 1.08] 0.016 1.001 1678 2163
Mean pressure utilized (partner’s view) 0.75* 0.09 [ 0.59, 0.96] 0.989 [0.92, 1.08] 0.044 1.002 1530 1753
Mean pushing experienced 0.60** 0.10 [ 0.43, 0.84] 0.998 [0.92, 1.08] 0.006 1.002 1717 2294
Mean pushing utilized (partner’s view) 0.74 0.13 [ 0.53, 1.03] 0.959 [0.92, 1.08] 0.082 1.002 1788 2445
Mean weartime NA NA NA NA NA NA NA NA NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced 0.99 0.04 [ 0.90, 1.08] 0.611 [0.84, 1.20] 1.000 1.000 5446 2698
Hu Daily persuasion utilized (partner’s view) 1.05 0.05 [ 0.96, 1.14] 0.845 [0.84, 1.20] 0.999 1.000 6950 2742
Hu Daily pressure experienced 1.19 0.13 [ 0.97, 1.49] 0.943 [0.84, 1.20] 0.508 1.002 6535 2884
Hu Daily pressure utilized (partner’s view) 1.11 0.12 [ 0.89, 1.38] 0.822 [0.84, 1.20] 0.745 1.002 7290 3194
Hu Daily pushing experienced 0.97 0.07 [ 0.84, 1.11] 0.665 [0.84, 1.20] 0.979 1.001 6965 3084
Hu Daily pushing utilized (partner’s view) 0.86* 0.06 [ 0.75, 0.99] 0.983 [0.84, 1.20] 0.669 1.000 7769 3417
Hu Day 1.20 0.16 [ 0.92, 1.56] 0.907 [0.84, 1.20] 0.500 1.002 7544 2807
Hu Daily weartime NA NA NA NA NA NA NA NA NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced 1.41 0.35 [ 0.87, 2.29] 0.916 [0.84, 1.20] 0.238 1.002 2418 2960
Hu Mean persuasion utilized (partner’s view) 1.48 0.36 [ 0.91, 2.41] 0.941 [0.84, 1.20] 0.195 1.002 2525 2587
Hu Mean pressure experienced 0.98 0.28 [ 0.55, 1.75] 0.525 [0.84, 1.20] 0.472 1.000 3786 2898
Hu Mean pressure utilized (partner’s view) 0.65 0.20 [ 0.36, 1.18] 0.920 [0.84, 1.20] 0.182 1.001 3543 2734
Hu Mean pushing experienced 0.57 0.22 [ 0.27, 1.22] 0.932 [0.84, 1.20] 0.136 1.002 2840 2878
Hu Mean pushing utilized (partner’s view) 0.72 0.28 [ 0.34, 1.53] 0.799 [0.84, 1.20] 0.256 1.002 2605 3077
Hu Mean weartime NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.30 0.04 [0.23, 0.41] NA NA NA 1.001 1167 1943
sd(Hurdle Intercept) 1.08 0.15 [0.83, 1.43] NA NA NA 1.001 981 1992
sd(Daily persuasion experienced) 0.11 0.02 [0.07, 0.16] NA NA NA 1.001 1649 2267
sd(Daily persuasion utilized (partner’s view)) 0.09 0.02 [0.05, 0.14] NA NA NA 1.005 1431 2110
sd(Daily pressure experienced) 0.06 0.06 [0.00, 0.22] NA NA NA 1.001 1164 2281
sd(Daily pressure utilized (partner’s view)) 0.05 0.04 [0.00, 0.18] NA NA NA 1.002 1286 2119
sd(Daily pushing experienced) 0.11 0.03 [0.06, 0.19] NA NA NA 1.002 1385 1569
sd(Daily pushing utilized (partner’s view)) 0.09 0.03 [0.03, 0.17] NA NA NA 1.002 1067 1295
sd(Hu Daily persuasion experienced) 0.48 0.08 [0.33, 0.67] NA NA NA 1.001 2363 2789
sd(Hu Daily persuasion utilized (partner’s view)) 0.38 0.08 [0.25, 0.57] NA NA NA 1.000 2347 3086
sd(Hu Daily pressure experienced) 0.29 0.24 [0.01, 0.91] NA NA NA 1.001 1282 2136
sd(Hu Daily pressure utilized (partner’s view)) 0.33 0.24 [0.02, 1.01] NA NA NA 1.001 1440 1566
sd(Hu Daily pushing experienced) 0.79 0.18 [0.50, 1.24] NA NA NA 1.002 2632 3185
sd(Hu Daily pushing utilized (partner’s view)) 0.77 0.18 [0.47, 1.20] NA NA NA 1.002 2671 2775
Additional Parameters
ar[1] 0.51 0.12 [0.34, 0.84] NA NA NA 1.192 15 28
sigma 0.35 0.06 [0.16, 0.45] NA NA NA 1.265 11 11 
# Plot continuous part of model

variable <- c(
  '(Intercept)',
  'b_persuasion_self_cw',
  'b_persuasion_partner_cw',
  'b_pressure_self_cw',
  'b_pressure_partner_cw',
  'b_pushing_self_cw',
  'b_pushing_partner_cw'
)


plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()

plot(
  bayestestR::rope(
    pa_sub, 
    parameter = variable, 
    range = rope_range_continuous,
    verbose = F,
    ci = 1
  )
) + theme_bw()

# Hurdle part of the model
variable <- c(
  'b_hu_persuasion_self_cw',
  'b_hu_persuasion_partner_cw',
  'b_hu_pressure_self_cw',
  'b_hu_pressure_partner_cw',
  'b_hu_pushing_self_cw',
  'b_hu_pushing_partner_cw'
)

plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

# The rope range for the bernoulli part of the model is -0.18, 0.18
plot(
  bayestestR::rope(pa_sub, parameter = variable, range = c(-0.18, 0.18), ci = 1),
  verbose = FALSE
) + theme_bw()

conds_eff <- conditional_spaghetti(
  pa_sub, 
  effects = c(
    'persuasion_self_cw',
    'persuasion_partner_cw',
    'pressure_self_cw',
    'pressure_partner_cw',
    'pushing_self_cw',
    'pushing_partner_cw'
  ),
  x_label = c(
    'Received Persuasion',
    'Exerted Persuasion',
    'Received Pressure',
    'Exerted Pressure',
    'Received Plan-Related Pushing',
    'Exerted Plan-Related Pushing'
  ),
  group_var = 'coupleID',
  plot_full_range = TRUE,
  y_limits = c(0, 100),
  y_label = "Same-Day MVPA",
  y_labels = c('Probability of Being Active', 'Minutes of MVPA When Active', 'Overall Expected Minutes of MVPA'),
  , filter_quantiles = .9995
  , font_family = 'Candara'
)
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
print(conds_eff)

$persuasion_self_cw

## Warning: Removed 99 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 73 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.00524

$persuasion_partner_cw

## Warning: Removed 105 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 78 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.00552

$pressure_self_cw

## Warning: Removed 58 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 53 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.015

$pressure_partner_cw

## Warning: Removed 15 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.0141

$pushing_self_cw

## Warning: Removed 58 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 37 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.00851

$pushing_partner_cw

## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.00807

Note. This graphic illustrates the relationship between social control and moderate to vigorous physical activity (MVPA) using a Bayesian Hurdle-Lognormal Multilevel Model. The predictor is centered within individuals to examine how deviations from their average social control relate to same-day MVPA. Shaded areas indicate credible intervals, thick lines show fixed effects, and thin lines represent random effects, highlighting variability across couples. The plots display the probability of being active, expected minutes of MVPA when active, and combined predicted MVPA. The bottom density plot visualizes the posterior distributions of slope estimates, transformed to represent multiplicative changes in odds ratios (hurdle component) or expected values. Medians and 95% credible intervals (2.5th and 97.5th percentiles) are shown. Effects are significant, when the 95% credible interval does not overlap 1.

Device Based MVPA

range(df_double$pa_obj, na.rm = T) 
## [1]   5.75 971.25
hist(df_double$pa_obj, breaks = 50)

df_double$pa_obj_log <- log(df_double$pa_obj)

hist(df_double$pa_obj_log, breaks = 50)

Lognormal Model

formula <- bf(
  pa_obj ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  , autocor = autocor_str
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 50)", class = "Intercept") 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
  , autocor_prior
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = lognormal()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = lognormal(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
pa_obj_log_digest <- digest::digest(pa_obj_log)
check_brms(pa_obj_log, log_pp_check = TRUE)
## 
## Divergences:
## 51 of 4000 iterations ended with a divergence (1.275%).
## Try increasing 'adapt_delta' to remove the divergences.
## 
## Tree depth:
## 1633 of 4000 iterations saturated the maximum tree depth of 10 (40.825%).
## Try increasing 'max_treedepth' to avoid saturation.
## 
## Energy:
## E-BFMI indicated possible pathological behavior:
##   Chain 1: E-BFMI = 0.029
##   Chain 2: E-BFMI = 0.012
##   Chain 3: E-BFMI = 0.010
##   Chain 4: E-BFMI = 0.055
## E-BFMI below 0.2 indicates you may need to reparameterize your model.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 2368 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## 
## Computed from 4000 by 3337 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -19618.6  90.2
## p_loo      3185.8  64.9
## looic     39237.2 180.4
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.0, 1.2]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)      969  29.0%   0       
##    (0.7, 1]   (bad)      1957  58.6%   <NA>    
##    (1, Inf)   (very bad)  411  12.3%   <NA>    
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(pa_obj_log, integer = TRUE, outliers_type = 'bootstrap')
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 245, observations = 3337, p-value <
## 2.2e-16
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.0005993407 0.0100614324
## sample estimates:
## outlier frequency (expected: 0.00383278393766856 ) 
##                                         0.07341924
# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)

summarize_brms(
  pa_obj_log, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  rope_range = rope_range_log,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
## Warning: There were 51 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning in summarize_brms(pa_obj_log, stats_to_report = c("CI", "SE", "pd", :
## Coefficients were exponentiated. Double check if this was intended.
exp(Est.) SE 95% CI pd ROPE inside ROPE Rhat Bulk_ESS Tail_ESS
Intercept 116.49*** 6.45 [104.27, 130.70] 1.000 [0.94, 1.07] 0.000 1.047 69 267
Within-Person Effects
Daily persuasion experienced 1.00 0.01 [ 0.98, 1.02] 0.555 [0.94, 1.07] 1.000 1.004 1073 719
Daily persuasion utilized (partner’s view) 1.00 0.01 [ 0.97, 1.02] 0.665 [0.94, 1.07] 1.000 1.001 1196 1929
Daily pressure experienced 1.04 0.03 [ 0.99, 1.10] 0.929 [0.94, 1.07] 0.818 1.002 1288 1874
Daily pressure utilized (partner’s view) 0.98 0.03 [ 0.93, 1.04] 0.719 [0.94, 1.07] 0.965 1.002 1841 2536
Daily pushing experienced 1.00 0.02 [ 0.97, 1.03] 0.552 [0.94, 1.07] 1.000 1.002 1702 2512
Daily pushing utilized (partner’s view) 1.02 0.02 [ 0.99, 1.06] 0.896 [0.94, 1.07] 0.995 1.003 1689 1619
Day 1.02 0.04 [ 0.95, 1.11] 0.731 [0.94, 1.07] 0.826 1.005 1609 2119
Daily weartime 1.00 0.00 [ 1.00, 1.00] 0.809 [0.94, 1.07] 1.000 1.001 2320 2657
Between-Person Effects
Mean persuasion experienced 1.05 0.08 [ 0.91, 1.21] 0.734 [0.94, 1.07] 0.515 1.007 673 1112
Mean persuasion utilized (partner’s view) 1.05 0.08 [ 0.90, 1.21] 0.726 [0.94, 1.07] 0.496 1.010 580 1397
Mean pressure experienced 0.85 0.07 [ 0.73, 1.00] 0.973 [0.94, 1.07] 0.116 1.003 1177 2076
Mean pressure utilized (partner’s view) 0.86 0.07 [ 0.74, 1.00] 0.974 [0.94, 1.07] 0.137 1.003 762 1745
Mean pushing experienced 1.14 0.12 [ 0.94, 1.39] 0.903 [0.94, 1.07] 0.226 1.008 960 1746
Mean pushing utilized (partner’s view) 1.14 0.11 [ 0.93, 1.37] 0.892 [0.94, 1.07] 0.244 1.009 788 1382
Mean weartime 1.00 0.00 [ 1.00, 1.00] 0.924 [0.94, 1.07] 1.000 1.001 1987 2612
Random Effects
sd(Intercept) 0.29 0.04 [0.23, 0.38] NA NA NA 1.006 609 1387
sd(Daily persuasion experienced) 0.06 0.01 [0.04, 0.09] NA NA NA 1.001 1102 1859
sd(Daily persuasion utilized (partner’s view)) 0.05 0.01 [0.03, 0.08] NA NA NA 1.003 974 1984
sd(Daily pressure experienced) 0.05 0.04 [0.00, 0.14] NA NA NA 1.005 580 1195
sd(Daily pressure utilized (partner’s view)) 0.03 0.03 [0.00, 0.11] NA NA NA 1.004 970 1715
sd(Daily pushing experienced) 0.07 0.03 [0.01, 0.14] NA NA NA 1.007 433 502
sd(Daily pushing utilized (partner’s view)) 0.03 0.02 [0.00, 0.09] NA NA NA 1.001 690 1564
Additional Parameters
ar[1] 0.39 0.11 [0.26, 0.71] NA NA NA 1.597 7 18
sigma 0.34 0.06 [0.20, 0.42] NA NA NA 1.647 7 20 
plot(
  bayestestR::p_direction(pa_obj_log),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(pa_obj_log, range = rope_range_log, ci = 1)
) + theme_bw()
## Possible multicollinearity between b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.76), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.7). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

# Nothing significant, no plots

Affect

range(df_double$aff, na.rm = T)
## [1] 0 5
hist(df_double$aff, breaks = 15)

Gaussian

formula <- bf(
  aff ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  , autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b")
  ,brms::set_prior("normal(0, 20)", class = "Intercept", lb=1, ub=6)
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
  , autocor_prior
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = gaussian()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood_gauss <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("mood_gauss", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
mood_gauss_digest <- digest::digest(mood_gauss)
check_brms(mood_gauss, log_pp_check = FALSE)
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## 
## Computed from 4000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -4820.4  63.5
## p_loo        94.1   4.4
## looic      9640.8 127.0
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.3]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(mood_gauss, integer = FALSE)

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 203, observations = 3736, p-value <
## 2.2e-16
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.04728469 0.06209462
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                           0.05433619
summarize_brms(
  mood_gauss, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F) %>%
  print_df(rows_to_pack = rows_to_pack)
Est. SE 95% CI pd ROPE inside ROPE Rhat Bulk_ESS Tail_ESS
Intercept 3.94*** 0.11 [ 3.72, 4.16] 1.000 [-0.11, 0.11] 0.000 1.006 642 1116
Within-Person Effects
Daily persuasion experienced 0.00 0.02 [-0.03, 0.03] 0.549 [-0.11, 0.11] 1.000 1.002 6740 3186
Daily persuasion utilized (partner’s view) -0.01 0.02 [-0.04, 0.02] 0.779 [-0.11, 0.11] 1.000 1.001 6310 3071
Daily pressure experienced 0.00 0.04 [-0.08, 0.07] 0.551 [-0.11, 0.11] 0.994 1.002 5819 3127
Daily pressure utilized (partner’s view) 0.05 0.04 [-0.03, 0.12] 0.874 [-0.11, 0.11] 0.960 1.001 6529 3445
Daily pushing experienced -0.03 0.02 [-0.08, 0.02] 0.913 [-0.11, 0.11] 1.000 1.002 6396 3213
Daily pushing utilized (partner’s view) -0.04 0.02 [-0.09, 0.01] 0.957 [-0.11, 0.11] 0.999 1.002 6118 3164
Day -0.09 0.08 [-0.26, 0.06] 0.882 [-0.11, 0.11] 0.593 1.002 7906 2488
Daily weartime NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 0.14 0.14 [-0.13, 0.41] 0.843 [-0.11, 0.11] 0.397 1.002 1435 2274
Mean persuasion utilized (partner’s view) 0.20 0.14 [-0.06, 0.48] 0.932 [-0.11, 0.11] 0.244 1.002 1445 2338
Mean pressure experienced -0.11 0.13 [-0.36, 0.13] 0.815 [-0.11, 0.11] 0.470 1.002 1835 2668
Mean pressure utilized (partner’s view) -0.01 0.13 [-0.28, 0.23] 0.545 [-0.11, 0.11] 0.636 1.001 2005 2496
Mean pushing experienced -0.19 0.20 [-0.58, 0.20] 0.831 [-0.11, 0.11] 0.288 1.000 2161 2661
Mean pushing utilized (partner’s view) -0.42* 0.20 [-0.81, -0.04] 0.984 [-0.11, 0.11] 0.057 1.001 2107 2757
Mean weartime NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.62 0.08 [0.49, 0.81] NA NA NA 1.003 896 1819
sd(Daily persuasion experienced) 0.03 0.02 [0.00, 0.07] NA NA NA 1.000 1930 2361
sd(Daily persuasion utilized (partner’s view)) 0.06 0.03 [0.01, 0.11] NA NA NA 1.006 991 1172
sd(Daily pressure experienced) 0.08 0.07 [0.00, 0.26] NA NA NA 1.001 1648 1773
sd(Daily pressure utilized (partner’s view)) 0.14 0.08 [0.01, 0.32] NA NA NA 1.003 1320 1550
sd(Daily pushing experienced) 0.09 0.03 [0.01, 0.16] NA NA NA 1.000 1562 895
sd(Daily pushing utilized (partner’s view)) 0.08 0.04 [0.01, 0.17] NA NA NA 1.003 1384 1721
Additional Parameters
ar[1] 0.45 0.02 [0.42, 0.48] NA NA NA 1.000 7873 3187
sigma 0.87 0.01 [0.85, 0.89] NA NA NA 1.000 7429 2816 
plot(
  bayestestR::p_direction(mood_gauss),
  priors = TRUE
)  + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(mood_gauss, ci = 1)
) + theme_bw()
## Possible multicollinearity between b_pressure_self_cb and
##   b_persuasion_self_cb (r = 0.8), b_pressure_partner_cb and
##   b_persuasion_self_cb (r = 0.75), b_pressure_self_cb and
##   b_persuasion_partner_cb (r = 0.76), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.79), b_pushing_partner_cb and
##   b_pushing_self_cb (r = 0.84). This might lead to inappropriate results.
##   See 'Details' in '?rope'.

conditional_spaghetti(
  mood_gauss, 
  effects = c('pushing_partner_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

$pushing_partner_cw

Reactance

range(df_double$reactance, na.rm = T) 
## [1] 0 5
hist(df_double$reactance, breaks = 7) 

hist(log(df_double$reactance+0.1), breaks = 10)

Ordinal

df_double$reactance_ordinal <- factor(df_double$reactance,
                                      levels = 0:5, 
                                      ordered = TRUE)

formula <- bf(
  reactance_ordinal ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  , autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , autocor_prior
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = cumulative() # HURDLE_CUMULATIVE
#  )


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance_ordinal <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::cumulative(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777
  , file = file.path("models_cache_brms", paste0("reactance_ordinal", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
reactance_ordinal_digest <- digest::digest(reactance_ordinal)
check_brms(reactance_ordinal)
## 
## Divergences:
## 21 of 4000 iterations ended with a divergence (0.525%).
## Try increasing 'adapt_delta' to remove the divergences.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 9 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## 
## Computed from 4000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -701.0 32.2
## p_loo       111.3  7.4
## looic      1402.0 64.3
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.1, 1.5]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     747   98.8%   46      
##    (0.7, 1]   (bad)        9    1.2%   <NA>    
##    (1, Inf)   (very bad)   0    0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(reactance_ordinal, outliers_type = 'bootstrap')
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 3, observations = 756, p-value = 0.02
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.002017196
## sample estimates:
## outlier frequency (expected: 0.000343915343915344 ) 
##                                         0.003968254
summarize_brms(
  reactance_ordinal, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_ordinal,
  model_rows_random = model_rows_random_ordinal,
  model_rownames_fixed = model_rownames_fixed_ordinal,
  model_rownames_random = model_rownames_random_ordinal,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack_ordinal)
## Warning: There were 21 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
OR SE 95% CI pd ROPE inside ROPE Rhat Bulk_ESS Tail_ESS
Intercepts
Intercept NA NA NA NA NA NA NA NA NA
Intercept[1] 2.47** 0.83 [ 1.29, 4.93] 0.997 [0.84, 1.20] 0.014 1.003 1520 1365
Intercept[2] 5.34*** 1.87 [ 2.71, 11.43] 1.000 [0.84, 1.20] 0.000 1.004 895 371
Intercept[3] 14.57*** 5.55 [ 7.25, 35.01] 1.000 [0.84, 1.20] 0.000 1.004 735 317
Intercept[4] 63.18*** 27.16 [ 29.03, 186.08] 1.000 [0.84, 1.20] 0.000 1.004 652 270
Intercept[5] 2363.89*** 1737.76 [673.98, 17307.95] 1.000 [0.84, 1.20] 0.000 1.004 615 294
Within-Person Effects
Daily persuasion experienced 1.04 0.07 [ 0.90, 1.19] 0.695 [0.84, 1.20] 0.982 1.002 4261 2229
Daily persuasion utilized (partner’s view) 0.92 0.07 [ 0.78, 1.08] 0.858 [0.84, 1.20] 0.871 1.001 5121 3041
Daily pressure experienced 0.88 0.11 [ 0.68, 1.12] 0.848 [0.84, 1.20] 0.645 1.002 3660 2026
Daily pressure utilized (partner’s view) 0.90 0.14 [ 0.65, 1.24] 0.748 [0.84, 1.20] 0.632 1.000 4002 1783
Daily pushing experienced 0.99 0.08 [ 0.84, 1.17] 0.537 [0.84, 1.20] 0.967 1.000 3991 2462
Daily pushing utilized (partner’s view) 1.07 0.11 [ 0.88, 1.32] 0.748 [0.84, 1.20] 0.847 1.000 3418 2059
Day 0.52 0.20 [ 0.24, 1.11] 0.956 [0.84, 1.20] 0.098 1.000 4260 3041
Daily weartime NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 0.89 0.39 [ 0.37, 2.15] 0.603 [0.84, 1.20] 0.310 1.001 2769 2658
Mean persuasion utilized (partner’s view) 1.12 0.55 [ 0.44, 2.96] 0.586 [0.84, 1.20] 0.286 1.002 2750 2252
Mean pressure experienced 1.94 0.90 [ 0.79, 4.82] 0.922 [0.84, 1.20] 0.112 1.001 3089 2730
Mean pressure utilized (partner’s view) 0.69 0.34 [ 0.26, 1.90] 0.771 [0.84, 1.20] 0.217 1.002 2456 1785
Mean pushing experienced 0.44 0.30 [ 0.11, 1.73] 0.887 [0.84, 1.20] 0.104 1.000 3745 3218
Mean pushing utilized (partner’s view) 1.15 0.81 [ 0.28, 4.62] 0.579 [0.84, 1.20] 0.195 1.000 3598 3423
Mean weartime NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 1.25 0.21 [ 0.92, 1.83] NA NA NA 1.002 1070 852
sd(Daily persuasion experienced) 0.20 0.14 [ 0.01, 0.50] NA NA NA 1.002 961 1569
sd(Daily persuasion utilized (partner’s view)) 0.16 0.13 [ 0.01, 0.49] NA NA NA 1.001 1304 1873
sd(Daily pressure experienced) 0.94 0.26 [ 0.54, 1.57] NA NA NA 1.001 1468 1393
sd(Daily pressure utilized (partner’s view)) 0.32 0.31 [ 0.02, 1.38] NA NA NA 1.002 916 1104
sd(Daily pushing experienced) 0.18 0.17 [ 0.01, 0.57] NA NA NA 1.003 865 1831
sd(Daily pushing utilized (partner’s view)) 0.11 0.11 [ 0.01, 0.53] NA NA NA 1.002 1588 2287
Additional Parameters
ar[1] 0.04 0.59 [-0.92, 0.93] NA NA NA 1.002 1533 2429
sigma NA NA NA NA NA NA NA NA NA
disc 1.00 0.00 [ 1.00, 1.00] NA NA NA NA NA NA 
plot(
  bayestestR::p_direction(reactance_ordinal),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()
## Warning: `b_Intercept[1]`, or one of its groups specified in `by`, is empty and
##   has no density information.
## Warning: `b_Intercept[2]`, or one of its groups specified in `by`, is empty and
##   has no density information.
## Warning: `b_Intercept[3]`, or one of its groups specified in `by`, is empty and
##   has no density information.
## Warning: `b_Intercept[4]`, or one of its groups specified in `by`, is empty and
##   has no density information.
## Warning: `b_Intercept[5]`, or one of its groups specified in `by`, is empty and
##   has no density information.

plot(
  bayestestR::rope(reactance_ordinal, range = c(-0.18, 0.18), ci = 1)
) + theme_bw()
## Possible multicollinearity between b_Intercept[4] and b_Intercept[2] (r
##   = 0.83), b_Intercept[4] and b_Intercept[3] (r = 0.9), b_Intercept[5] and
##   b_Intercept[4] (r = 0.7), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.72). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

conditional_spaghetti(
  reactance_ordinal, 
  effects = c('persuasion_self_cw', 'pressure_self_cw')
  , group_var = 'coupleID'
  #, n_groups = 15
  , plot_full_range = T
)

\(persuasion_self_cw <img src="02_SensitivityAr1_files/figure-html/report_reactance_ordinal-3.png" width="2400" />\)pressure_self_cw

Binary

introduce_binary_reactance <- function(data) {
  data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
  return(data)
}



df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
  for (i in seq_along(implist)) {
    implist[[i]] <- introduce_binary_reactance(implist[[i]])
  }
}


formula <- bf(
  is_reactance ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  , autocor = autocor_str
  )



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5) 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , autocor_prior
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = bernoulli()
#  )



#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::bernoulli(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("is_reactance", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
is_reactance_digest <- digest::digest(is_reactance)
check_brms(is_reactance)
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 744 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## 
## Computed from 4000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -361.7 13.2
## p_loo       299.6 11.5
## looic       723.3 26.4
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 1.3]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)      12    1.6%   300     
##    (0.7, 1]   (bad)      387   51.2%   <NA>    
##    (1, Inf)   (very bad) 357   47.2%   <NA>    
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(is_reactance, integer = FALSE)

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 7, observations = 756, p-value = 0.0009473
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.003730576 0.018984066
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.009259259
summarize_brms(
  is_reactance,
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
OR SE 95% CI pd ROPE inside ROPE Rhat Bulk_ESS Tail_ESS
Intercept 2.15 1.85 [0.40, 13.28] 0.819 [0.83, 1.20] 0.122 1.001 2048 2467
Within-Person Effects
Daily persuasion experienced 1.15 0.36 [0.64, 2.20] 0.684 [0.83, 1.20] 0.412 1.001 1332 2080
Daily persuasion utilized (partner’s view) 0.97 0.35 [0.44, 2.16] 0.531 [0.83, 1.20] 0.381 1.000 1604 2069
Daily pressure experienced 0.84 0.45 [0.28, 2.50] 0.621 [0.83, 1.20] 0.255 1.003 2016 2080
Daily pressure utilized (partner’s view) 0.94 0.69 [0.20, 4.75] 0.531 [0.83, 1.20] 0.190 1.002 1641 1946
Daily pushing experienced 0.98 0.34 [0.46, 1.96] 0.524 [0.83, 1.20] 0.401 1.001 1761 2271
Daily pushing utilized (partner’s view) 0.97 0.47 [0.38, 2.68] 0.521 [0.83, 1.20] 0.298 1.002 1924 2119
Day 0.19 0.26 [0.01, 3.18] 0.887 [0.83, 1.20] 0.047 1.002 2111 2578
Daily weartime NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 1.09 1.54 [0.07, 14.63] 0.525 [0.83, 1.20] 0.108 1.001 1686 2507
Mean persuasion utilized (partner’s view) 1.48 2.08 [0.10, 28.28] 0.616 [0.83, 1.20] 0.103 1.001 2205 2771
Mean pressure experienced 1.81 2.68 [0.10, 36.78] 0.665 [0.83, 1.20] 0.099 1.004 1857 2516
Mean pressure utilized (partner’s view) 0.88 1.29 [0.05, 15.37] 0.536 [0.83, 1.20] 0.104 1.001 2194 2614
Mean pushing experienced 0.13 0.25 [0.00, 6.05] 0.860 [0.83, 1.20] 0.046 1.002 1996 2690
Mean pushing utilized (partner’s view) 2.43 4.74 [0.05, 118.99] 0.682 [0.83, 1.20] 0.064 1.001 2543 2902
Mean weartime NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 5.94 1.19 [ 3.83, 8.54] NA NA NA 1.002 1092 1437
sd(Daily persuasion experienced) 1.15 0.59 [ 0.16, 2.55] NA NA NA 1.010 405 717
sd(Daily persuasion utilized (partner’s view)) 1.08 0.70 [ 0.08, 2.74] NA NA NA 1.016 517 898
sd(Daily pressure experienced) 3.11 0.97 [ 1.50, 5.26] NA NA NA 1.004 1196 1439
sd(Daily pressure utilized (partner’s view)) 0.97 0.87 [ 0.04, 3.64] NA NA NA 1.001 1241 1641
sd(Daily pushing experienced) 0.71 0.59 [ 0.03, 2.11] NA NA NA 1.001 685 1431
sd(Daily pushing utilized (partner’s view)) 0.62 0.54 [ 0.03, 2.24] NA NA NA 1.003 1164 2081
Additional Parameters
ar[1] 0.11 0.10 [-0.09, 0.32] NA NA NA 1.011 584 1021
sigma NA NA NA NA NA NA NA NA NA 
plot(
  bayestestR::p_direction(is_reactance),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(is_reactance, ci = 1)
) + theme_bw()

conditional_spaghetti(
  is_reactance, 
  effects = c('pressure_self_cw', 'pushing_self_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

\(pressure_self_cw <img src="02_SensitivityAr1_files/figure-html/report_is_reactance-3.png" width="2400" />\)pushing_self_cw

Report All Models

process_model_component <- function(obj) {
  # Convert to string, modify, and evaluate
  name <- deparse(substitute(obj))
  if (report_hurdle) name <- paste0(name, '_hu')
  if (report_ordinal) name <- paste0(name, '_ordinal')
  return(get(name, envir = parent.frame()))
}



all_models <- report_side_by_side(
  pa_sub,
  pa_obj_log,
  mood_gauss,
  reactance_ordinal,
  is_reactance,
  
  stats_to_report = c('CI', 'pd'),
  
  model_rows_random = process_model_component(model_rows_random),
  model_rows_fixed = process_model_component(model_rows_fixed),
  model_rownames_random = process_model_component(model_rownames_random),
  model_rownames_fixed = process_model_component(model_rownames_fixed)
)

[1] “pa_sub”

## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “pa_obj_log”

## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
## Warning: There were 51 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “mood_gauss” [1] “reactance_ordinal”

## Warning: There were 21 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

[1] “is_reactance”

# pretty printing

summary_all_models <- all_models %>%
  print_df(rows_to_pack = process_model_component(rows_to_pack)) %>%
  add_header_above(
    c(" ", "Subjective MVPA Hurdle Lognormal" = (length(all_models) / 5),  
      "Device-Based MVPA Log (Gaussian)" = (length(all_models) / 5), 
      "Mood Gaussian" = (length(all_models) / 5),
      "Reactance Ordinal" = (length(all_models) / 5),
      "Reactance Dichotome" = (length(all_models) / 5))
  )


export_xlsx(
  summary_all_models, 
  rows_to_pack = process_model_component(rows_to_pack),
  file.path("Output", "AllModelsWithAR1.xlsx"), 
  merge_option = 'header', 
  simplify_2nd_row = TRUE,
  line_above_rows = c(1,2),
  line_below_rows = c(-1)
)

  
print(summary_all_models)
Subjective MVPA Hurdle Lognormal
Device-Based MVPA Log (Gaussian)
Mood Gaussian
Reactance Ordinal
Reactance Dichotome
exp(Est.) pa_sub 95% CI pa_sub pd pa_sub exp(Est.) pa_obj_log 95% CI pa_obj_log pd pa_obj_log Est. mood_gauss 95% CI mood_gauss pd mood_gauss OR reactance_ordinal 95% CI reactance_ordinal pd reactance_ordinal OR is_reactance 95% CI is_reactance pd is_reactance
Intercept 44.10*** [38.98, 49.78] 1.000 116.49*** [104.27, 130.70] 1.000 3.94*** [ 3.72, 4.16] 1.000 NA NA NA 2.15 [0.40, 13.28] 0.819
Hurdle Intercept 0.89 [ 0.51, 1.67] 0.654 NA NA NA NA NA NA NA NA NA NA NA NA
Conditional Within-Person Effects
Daily persuasion experienced 0.97 [ 0.93, 1.01] 0.938 1.00 [ 0.98, 1.02] 0.555 0.00 [-0.03, 0.03] 0.549 1.04 [ 0.90, 1.19] 0.695 1.15 [0.64, 2.20] 0.684
Daily persuasion utilized (partner’s view) 1.00 [ 0.96, 1.04] 0.504 1.00 [ 0.97, 1.02] 0.665 -0.01 [-0.04, 0.02] 0.779 0.92 [ 0.78, 1.08] 0.858 0.97 [0.44, 2.16] 0.531
Daily pressure experienced 1.06 [ 0.96, 1.16] 0.882 1.04 [ 0.99, 1.10] 0.929 0.00 [-0.08, 0.07] 0.551 0.88 [ 0.68, 1.12] 0.848 0.84 [0.28, 2.50] 0.621
Daily pressure utilized (partner’s view) 0.99 [ 0.90, 1.08] 0.615 0.98 [ 0.93, 1.04] 0.719 0.05 [-0.03, 0.12] 0.874 0.90 [ 0.65, 1.24] 0.748 0.94 [0.20, 4.75] 0.531
Daily pushing experienced 1.03 [ 0.97, 1.09] 0.845 1.00 [ 0.97, 1.03] 0.552 -0.03 [-0.08, 0.02] 0.913 0.99 [ 0.84, 1.17] 0.537 0.98 [0.46, 1.96] 0.524
Daily pushing utilized (partner’s view) 0.99 [ 0.93, 1.05] 0.657 1.02 [ 0.99, 1.06] 0.896 -0.04 [-0.09, 0.01] 0.957 1.07 [ 0.88, 1.32] 0.748 0.97 [0.38, 2.68] 0.521
Day 1.13 [ 1.00, 1.29] 0.972 1.02 [ 0.95, 1.11] 0.731 -0.09 [-0.26, 0.06] 0.882 0.52 [ 0.24, 1.11] 0.956 0.19 [0.01, 3.18] 0.887
Daily weartime NA NA NA 1.00 [ 1.00, 1.00] 0.809 NA NA NA NA NA NA NA NA NA
Conditional Between-Person Effects
Mean persuasion experienced 1.07 [ 0.87, 1.31] 0.746 1.05 [ 0.91, 1.21] 0.734 0.14 [-0.13, 0.41] 0.843 0.89 [ 0.37, 2.15] 0.603 1.09 [0.07, 14.63] 0.525
Mean persuasion utilized (partner’s view) 1.14 [ 0.94, 1.40] 0.908 1.05 [ 0.90, 1.21] 0.726 0.20 [-0.06, 0.48] 0.932 1.12 [ 0.44, 2.96] 0.586 1.48 [0.10, 28.28] 0.616
Mean pressure experienced 1.39** [ 1.10, 1.79] 0.997 0.85 [ 0.73, 1.00] 0.973 -0.11 [-0.36, 0.13] 0.815 1.94 [ 0.79, 4.82] 0.922 1.81 [0.10, 36.78] 0.665
Mean pressure utilized (partner’s view) 0.75* [ 0.59, 0.96] 0.989 0.86 [ 0.74, 1.00] 0.974 -0.01 [-0.28, 0.23] 0.545 0.69 [ 0.26, 1.90] 0.771 0.88 [0.05, 15.37] 0.536
Mean pushing experienced 0.60** [ 0.43, 0.84] 0.998 1.14 [ 0.94, 1.39] 0.903 -0.19 [-0.58, 0.20] 0.831 0.44 [ 0.11, 1.73] 0.887 0.13 [0.00, 6.05] 0.860
Mean pushing utilized (partner’s view) 0.74 [ 0.53, 1.03] 0.959 1.14 [ 0.93, 1.37] 0.892 -0.42* [-0.81, -0.04] 0.984 1.15 [ 0.28, 4.62] 0.579 2.43 [0.05, 118.99] 0.682
Mean weartime NA NA NA 1.00 [ 1.00, 1.00] 0.924 NA NA NA NA NA NA NA NA NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced 0.99 [ 0.90, 1.08] 0.611 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily persuasion utilized (partner’s view) 1.05 [ 0.96, 1.14] 0.845 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pressure experienced 1.19 [ 0.97, 1.49] 0.943 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pressure utilized (partner’s view) 1.11 [ 0.89, 1.38] 0.822 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pushing experienced 0.97 [ 0.84, 1.11] 0.665 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pushing utilized (partner’s view) 0.86* [ 0.75, 0.99] 0.983 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Day 1.20 [ 0.92, 1.56] 0.907 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily weartime NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced 1.41 [ 0.87, 2.29] 0.916 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean persuasion utilized (partner’s view) 1.48 [ 0.91, 2.41] 0.941 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pressure experienced 0.98 [ 0.55, 1.75] 0.525 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pressure utilized (partner’s view) 0.65 [ 0.36, 1.18] 0.920 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pushing experienced 0.57 [ 0.27, 1.22] 0.932 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pushing utilized (partner’s view) 0.72 [ 0.34, 1.53] 0.799 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean weartime NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.30 [0.23, 0.41] NA 0.29 [0.23, 0.38] NA 0.62 [0.49, 0.81] NA 1.25 [ 0.92, 1.83] NA 5.94 [ 3.83, 8.54] NA
sd(Hurdle Intercept) 1.08 [0.83, 1.43] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Daily persuasion experienced) 0.11 [0.07, 0.16] NA 0.06 [0.04, 0.09] NA 0.03 [0.00, 0.07] NA 0.20 [ 0.01, 0.50] NA 1.15 [ 0.16, 2.55] NA
sd(Daily persuasion utilized (partner’s view)) 0.09 [0.05, 0.14] NA 0.05 [0.03, 0.08] NA 0.06 [0.01, 0.11] NA 0.16 [ 0.01, 0.49] NA 1.08 [ 0.08, 2.74] NA
sd(Daily pressure experienced) 0.06 [0.00, 0.22] NA 0.05 [0.00, 0.14] NA 0.08 [0.00, 0.26] NA 0.94 [ 0.54, 1.57] NA 3.11 [ 1.50, 5.26] NA
sd(Daily pressure utilized (partner’s view)) 0.05 [0.00, 0.18] NA 0.03 [0.00, 0.11] NA 0.14 [0.01, 0.32] NA 0.32 [ 0.02, 1.38] NA 0.97 [ 0.04, 3.64] NA
sd(Daily pushing experienced) 0.11 [0.06, 0.19] NA 0.07 [0.01, 0.14] NA 0.09 [0.01, 0.16] NA 0.18 [ 0.01, 0.57] NA 0.71 [ 0.03, 2.11] NA
sd(Daily pushing utilized (partner’s view)) 0.09 [0.03, 0.17] NA 0.03 [0.00, 0.09] NA 0.08 [0.01, 0.17] NA 0.11 [ 0.01, 0.53] NA 0.62 [ 0.03, 2.24] NA
sd(Hu Daily persuasion experienced) 0.48 [0.33, 0.67] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily persuasion utilized (partner’s view)) 0.38 [0.25, 0.57] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pressure experienced) 0.29 [0.01, 0.91] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pressure utilized (partner’s view)) 0.33 [0.02, 1.01] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pushing experienced) 0.79 [0.50, 1.24] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pushing utilized (partner’s view)) 0.77 [0.47, 1.20] NA NA NA NA NA NA NA NA NA NA NA NA NA
Additional Parameters
ar[1] 0.51 [0.34, 0.84] NA 0.39 [0.26, 0.71] NA 0.45 [0.42, 0.48] NA 0.04 [-0.92, 0.93] NA 0.11 [-0.09, 0.32] NA
sigma 0.35 [0.16, 0.45] NA 0.34 [0.20, 0.42] NA 0.87 [0.85, 0.89] NA NA NA NA NA NA NA 
report::report_system()

Analyses were conducted using the R Statistical language (version 4.4.1; R Core Team, 2024) on Windows 11 x64 (build 22635)

report::report_packages()
  • rmarkdown (version 2.29; Allaire J et al., 2024)
  • beepr (version 2.0; Bååth R, 2024)
  • R.methodsS3 (version 1.8.2; Bengtsson H, 2003)
  • R.oo (version 1.27.0; Bengtsson H, 2003)
  • R.utils (version 2.12.3; Bengtsson H, 2023)
  • brms (version 2.22.0; Bürkner P, 2017)
  • posterior (version 1.6.0; Bürkner P et al., 2024)
  • extrafont (version 0.19; Chang W, 2023)
  • digest (version 0.6.37; Eddelbuettel D, 2024)
  • Rcpp (version 1.0.13.1; Eddelbuettel D et al., 2024)
  • bayesplot (version 1.11.1; Gabry J, Mahr T, 2024)
  • lubridate (version 1.9.3; Grolemund G, Wickham H, 2011)
  • DHARMa (version 0.4.7; Hartig F, 2024)
  • fs (version 1.6.5; Hester J et al., 2024)
  • wbCorr (version 0.1.22; Küng P, 2023)
  • see (version 0.9.0; Lüdecke D et al., 2021)
  • tibble (version 3.2.1; Müller K, Wickham H, 2023)
  • patchwork (version 1.3.0; Pedersen T, 2024)
  • R (version 4.4.1; R Core Team, 2024)
  • openxlsx (version 4.2.7.1; Schauberger P, Walker A, 2024)
  • ggplot2 (version 3.5.1; Wickham H, 2016)
  • forcats (version 1.0.0; Wickham H, 2023)
  • stringr (version 1.5.1; Wickham H, 2023)
  • rvest (version 1.0.4; Wickham H, 2024)
  • tidyverse (version 2.0.0; Wickham H et al., 2019)
  • readxl (version 1.4.3; Wickham H, Bryan J, 2023)
  • dplyr (version 1.1.4; Wickham H et al., 2023)
  • purrr (version 1.0.2; Wickham H, Henry L, 2023)
  • readr (version 2.1.5; Wickham H et al., 2024)
  • xml2 (version 1.3.6; Wickham H et al., 2023)
  • tidyr (version 1.3.1; Wickham H et al., 2024)
  • ggridges (version 0.5.6; Wilke C, 2024)
  • knitr (version 1.49; Xie Y, 2024)
  • kableExtra (version 1.4.0; Zhu H, 2024)
report::cite_packages()